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Abstract 

This study presents a transient nonlinear finite element analysis within the realm of 
a multi-body dynamics formulation for determining the dynamic response of a 
moderately thick laminated shell undergoing a rapid and large rotational motion and 
nonlinear elastic deformations. Nonlinear strain measure and rotation, as well as the 
transverse shear deformation, are explicitly included in the formulation in order to 
capture the proper motion-induced stiffness of the laminate. The equations of motion 
are derived from the virtual work principle. The analysis utilizes a shear deformable 
shallow shell element along with the co-rotational form of the updated Lagrangian 
formulation. The shallow shell element formulation is based on the Reissner-Mindlin and 
Marguerre theory. 
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Introduction 

Thick laminated composites such as helicopter blades can experience large elastic 
displacement and rotations due to rapid and large rotational motions. Any attempt to 
design composite helicopter blades must include the effect of motion-induced stiffness 
in both the longitudinal and chordwise directions. The stress analysis capability must 
account for the interaction of the nonlinear elastic deformations with the overall 
dynamic motion. 

Review of previous studies relevant to helicopter blade design reveals that the 
majority of the analyses are limited to beam-type approximations, disregarding the 
chordwise deformations. In order to eliminate this shortcoming, Bhumbla and Kosmatka 
[ 1 ]. Bauchau and Chiang [2], and Kosmatka [3] employed shear deformable shell 
elements within analyses restricted to rotating laminates under constant angular 
velocity. Their formulations exclude the effect of spin-up time, during which the angular 
velocity is not constant. Furthermore, the influence of the coupling between the rigid- 
body motion and the flexible deformation is not reflected in the construction of the 
mass matrix in any of these models. The inertial loads are treated as known external 
forces. Within multibody dynamics formalism, the mass matrix should include the 
coupling between the rigid-body motion, inertial forces, and the elastic deformations so 
that the inertial loads can be determined as part of the solution. Changes in inertial 
loads due to spin-up or external loads become important, especially in aeroelastic 
tailoring of the composite blades. 

In order to account for this discrepancy, Bauchau and Kang [4] considered a 
multibody dynamics formulation coupled with elastic deformations. Their analysis, 
however, employs only beam-type finite elements. Extensive reviews of relevant 
multibody dynamics formulations coupled with elastic deformations can be found in 
Boutaghou et al. [5] and Tsang [6], 

Although the response of an elastic beam subjected to a sudden slew rate coupled 
with large elastic deformations has been investigated extensively, there are no analyses 
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regarding the dynamic response of an elastic shell undergoing large overall motion 
accompanied by large elastic displacements and rotations. Banarjee and Kane [7], 
Banerjee and Dickens [8], Chang and Shabana [ 9—1 0] , Boutaghou et al. [ 5] , and Tsang 
[6] included the interaction of large overall motion and small strain elastic 
deformations. The complexity of the problem arising from the interaction of large rigid- 
body translational and rotational motions and large displacements and rotations renders 
the analytical and numerical predictions very difficult. 

All of these previous analyses are limited to isotropic plates and suffer from loss 
of additional motion-induced stiffness due to large elastic displacements (membrane 
and bending) and rotations. Also, they require that the rigid-body motion be known a 
priori . The coupling between bending and stretching and bending and twisting in 
anisotropic composite plates and the presence of initial curvature further complicate the 
problem. 

The present analysis is aimed at capturing the motion-induced membrane and 
bending stiffness variations by explicitly including the nonlinear strain measure and 
rotations and the transverse shear deformations. It is to be noted that none of the 
available finite element analyses, including commercial programs, have the capability to 
solve this problem within the realm of multibody dynamics formalism. The transverse 
shear deformation, which is essential for accurate failure prediction of thick composite 
laminates, is included by extending Tessler's [ 1 1 ] shallow shell element into the 
nonlinear range and its implementation into the transient finite element formulation 
within the framework of multibody dynamics formalism. 

Equations of Motion 

The kinematic description of a point P in a three-noded shallow shell element 
between time t and t + At is illustrated in Fig. 1, The equations of motion are derived 
based on the concept of virtual displacements in conjunction with the co-rotational form 
of the updated Lagrangian description of motion. The formulation utilizes the inertial, 
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body-fixed, and element Cartesian coordinates as (x, y, z), (x' t y', z'), and (x", y", z M ), 
respectively. The configuration at time t + At represents the unknown current 
equilibrium configuration. The last known equilibrium configuration is at time t, and it is 
referred to as the updated configuration. In the updated Lagrangian form of the co- 
rotational approach, the current configuration utilizes the projection of the deformed 
element surface onto the x M -y" plane of the element coordinate system in the updated 
configuration as a reference configuration. The angular velocity, Q(t), of the shell 
arising from either a prescribed slew rate or external forces results in the rigid-body 
motion of the element coordinate system as the deformation proceeds. 

Matrices Z and T are utilized in transforming a vector from the element coordinate 
system to the body-fixed frame and from the body-fixed frame to the inertial frame, 
respectively. A vector with a single prime or a double prime is defined with respect to 
base vectors of the body-fixed or element coordinate system, respectively. 

In the absence of body forces, the virtual work done at time t + At in the kth 
element due to inertial, internal, and external forces is expressed as 

5 t+At R T UAt R d l v + 

Jv 

In this expression, the left subscript indicates the configuration by which the quantity is 
measured. The left superscript refers to the configuration of the body at a specific 
time. The volume of the kth facet shell element containing the point P at time t is 
denoted by V. Differentiation with respect to time is denoted by (•) superscribed above 
the variable, and its virtual value by 8, 

The position of point P at time t -h At is described by the vector u ^ l R in reference 
to the origin of the inertial coordinate system, (x, y z). The components of the 
incremental Green strain tensor and the Pioia-Kirchhoff stress tensor associated with the 
kth element are contained in vectors t E and respectively. The virtual work for the 

kth element due to nodal forces and moments at time t + At is denoted by 8W . The 
constant mass density of the material is denoted by p. 


S t E' 


t+A !s d l V 


8W 


( 1 ) 
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Virtual Work by Inertial Forces. Based on the kinematics of motion as 
illustrated in Fig. 1, the position of point P in the kth element at time t + At is 
expressed as 

UAt R . t + At R() + U4t r ( 2) 

where ,+At R 0 describes the rigid-body motion of the shell by the position of the origin 
of the body-fixed frame (x', y', z') at time t + At. Relative to this body-fixed frame, the 
point P undergoes rigid- and flexible-body motions from the previously known 
configuration to the current (unknown) configuration. Between time t and t + At, the 
rigid-body motion and elastic deformation of point P are specified by l r and t+At u, 
respectively. In reference to the body-fixed frame, the position of point P at time t + At 
is given by t+At r. (In the remaining expressions, the left superscript is dropped for 
simplicity unless it is different from t + At.) Applying the variational operator to Eq. 
(2), with the use of the transformation 

r - Tr’ , with r' - l r' + u' (3) 

and the variation of T given by 

5T = 5^ T (4) 

puts Eq. (2) in the following form: 

5R - 5R 0 - T V T T hw- + T6u' (5) 

where &w is the virtual rigid-body rotation vector. Hereafter, indicates a skew- 
symmetric matrix and hw- is a skew-symmetric matrix with virtual rotations of the body- 
fixed frame. 

The acceleration at point P is given by 


R - R 0 - TT T T Q + Tu' + 2QTu‘ + QQTr' 


( 6 ) 
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Relative to the body-fixed frame, the incremental displacement vector, u', is 
expressed in terms of the incremental total nodal vector, v', as 


u’ - N>’ 


(7) 


based on the finite element displacement approximation with the interpolation matrix, 
N^. The total nodal vector consists of the incremental elastic displacement, u' 0i , and 
rotation, 0' Oi , at each node of the kth shallow shell element shown in Fig. 2 and is 
defined as 

v ,T - {v? v'J Vj) (8) 


in which v'J - {u’J 0'J}. As suggested by Christensen and Lee [ 12], the negligible time 
variation of the interpolation matrix, NJ^, leads to 


u’ - N^v' and ii - N' m v' 


(9) 


Utilizing the transformation matrix, Z, associated with the kth element, the incremental 
displacement vector, u', is expressed as 


u' - Zu" (10) 

The transformation between the incremental total nodal vectors defined relative to the 
body-fixed and element frames is achieved by 

v’ - Z'v" (11) 

The transformation matrix, Z‘, is defined as 

*ii iJ - 1.3 (12) 

in which 5y is the Kronecker delta and Z 9 is the transformation matrix relating the 
rotations between the element coordinate system and the body-fixed frame as 
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“ ZflG'b (13) 

Using the interpolation matrix, N'^, the incremental displacement vector, u", is 
approximated in terms of the incremental total nodal vector, v", as 


u 


ii 


N" y" 
™ m ¥ 


(14) 


Based on the Reissner-Mindlin plate theory, the incremental displacement vector 
for a flat laminate can be expressed as 



r 

*(x'\ y", z") 


-o(x". y") 


V( x "* y") 

u M - h 

*(x", y". z") 

► » < 

-o(x". y") 

- + z M < 

®x"( x ". y") > 


Wx". y". z") 


-o(x". y") _ 


0 

t J 


(15a) 


or 


u 


* i 


+ z 


*®o 


(15b) 


in which a 0t ^ 0 , and represent mid-surface displacements. The in-plane 

displacement components in the x" and y" directions are denoted by ^ and 
respectively. The out-of-plane deflection is represented by a*. The incremental 
displacement and rotation vectors associated with the mid-surface, u" 0 and fl'b, are 
approximated in terms of the nodal unknown displacements as 


“ o = ^mu v " and 0” 0 - N^v" (16) 

where v" is the incremental total nodal vector and and N m $ are the interpolation 
matrices for translational and rotational displacement fields. Substituting for the mid- 
surface displacement and rotation vectors, u’q and ]n Eq. (15) results in 

u " * (^mu + Z M N mfl )v M (17) 

Thus, the explicit form of the interpolation matrix, N'^, in Eq. (14) can be established 


as 
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N; “ Nmo + z"Nmfl 


(18) 


Substituting from Eq. (14) into (10) in conjunction with Eq. (11) results in 


u' - ZN'^Z* 7 v' 


(19) 


Comparison of Eqs. (19) and (7) provides the expression for NL as 


Nm - ZN’InZ- 1 


( 20 ) 


After substituting from Eq. (9), along with Eq. (20), into Eqs. (5) and (6), the first 
integral, corresponding to the virtual work due to inertial forces in the kth element at 
time t + At, can be expressed in matrix form as 


1 . 


p 5R T RdV - 
V 


5R 0 

bw t 

by' 


M RR M Ra , M Rv . 
M 

sym 




in which 

Mrr = P 'V I , 


M vV 


-] 

* ^ 


* N 

T 

r 

1 


Ro 


5R 0 



+ \ 2 

< 

Q 

► + « 

blD- 

► < 

9i 

+ 92 


v’ 


5v' 


h i 

+ h 2 


(21) 


m r^ " -PT 


r'dV T t 


m Rv' ” pTZ 


J *V 


N'JpdV Z* T 


~ -pTZ 


r” r"dV Z T T T 


M^ v . - pTZ 


r-N'^dV Z* T 


M vV - pZ [ N" m J N" m d l W Z‘ T 
J l V 


f, - 2pTZ£2" 


N'mdV V" , 


f 2 - pTQ'Q’ r'dV 

J*V 


9i - 2pTZ 


r M Q M N'^d , V v" , 


g 2 - pTZ 


r*'fi ,, Q ,, r"d t V 


V 
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hi 


2pZ* 



„ T 
m 


Q"N 


it 

m 


d l V v", 


oz - [ N'^ T 
J'v 


Q"Q"r"d V 


where the identity matrix is denoted by I. 

Virtual work by internal stresses. The second integral, representing the virtual 
work by internal stresses, in Eq. (1) can be expressed for a laminated composite 
material in the form [ 1 3] 



«tEja 


rfr) 


A 


dz d l a + 



V k > 


dz d l a 


+ X ft f 5 *El t ^ (k> dz d l a (22) 

M J A Vi 


where K denotes the number of layers in the laminate. The position of the layers in 
reference to the mid-surface is specified by h k . The thickness of the kth layer is given 
by t k - h k - h k-1 . At time t, l A represents the projections of the surface area on the 
(x'\ y") plane. The stress and strain components are related through the relation a * 

E L , in which is for the kth orthotropic lamina referenced to the arbitrary axes. 
Each layer is assumed to be homogeneous, elastic, and orthotropic with elastic moduli, 
E l and E t ; shear modulus, G LT ; and Poisson's ratio, v LJ . The subscripts L and T 
specify the longitudinal and transverse directions relative to the fibers in the layer with 
thickness t. 

The linear and nonlinear components of the Green's strain tensor in vector form are 
expressed as 
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u '0,x n 


e y”.x" 


' 

h 0 x - 

4 ^0,y" 

' + z" - 

0 X „ y» 

► _ < 

^0,y**^x" * 

^0,y H+ ^0,x" 


0 


^0,x"®x" ^ ^O.y^y" 


Enl- 



2 2 2 

“'O.x" + ^O.x” + “'‘0.x" 

2 2 2 

0-O.y" + ^0,/* + “*0,y" 


2(^o,x“^O.y" + ^O.x’^O.y" + “'"0. x'^O.y") 


( 23 ) 


where h 0 (x", y") describes the surface of the shallow shell. As given by Tessler ( 1 1 ] , 
the linear component of the Green strain tensor is based on the Reissner-Mindlin 
definition combined with the Marguerre shallow shell theory. The nonlinear strain 
vector, E^, is an approximation to its complete form. 

Substituting for the strain components and performing integration along the 
transverse direction in Eq. (22) lead to 


■ 

■ r 

5 t e b C t e L d'a + 

J'a J 

. S t E ^L * s m d * a + 

'a J 


6 t eL *s 


d'a 


(24) 


in which C is composed of the extensional, coupling, bending, and transverse stiffness 
matrices, A, B, D, and G 


C 


A 

B 

0 


B 0 

D 0 

0 G 


the vector l s, which contains the in-plane resultant forces, t s m , moments, *s b , and shear 
resultant forces, l s q , is in the form 


tsT = i' s m 


ts b ts ql 


(25) 
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The vector t e L is composed of incremental mid-plane strain, curvature, and transverse 
shear components, 


t®L - {a-Q,*" “ ^0,y” ” ^O.y“^x"' “"O.y” + ^O.x" ” ^0,x“®x" 

— h O y -0y., ^y" ^", 0 X " y» , ® x " x" ^y“.y "• “^o.x" ®y"» ^o.y” ®x"l (26) 


The incremental in-plane mid-plane displacement components are approximated by 
cubic interpolation functions, and -^' ac , as in Tessler ( 1 1 ], 


r ^ 

USQ 

< 

— V 

yc 

A 

* ^ 
^0k 

9 

,+ 

"Ok 

"0 

V. J 

k-1 

_ ^0k 

k-4 



(27) 


where u-^ and v ^ are the in-plane displacements specified at the nodes. The 
displacements corresponding to the corner nodes are specified by a 0k , ^ 0k with k - 1, 
2, 3. The remaining » 0k and ^ represent the intra-edge displacements. The 
displacements at the center of the element are given by u^. The mid-plane transverse 
displacement field is approximated by a quadratic interpolation function, Jf sk , as 


o 


6 

y. -^sk^Ok 
k-1 


(28) 


where correspond to the corner and mid-edge nodes. 

The mid-plane out-of-plane rotational components are approximated by linear 
interpolation functions, £ k> as 



(29) 


where 0 X .^ and 0 y » k are nodal rotations at the corner nodes. 
In matrix notation, Eqs. (27)— (29) can be rewritten as 


u 


* i 
0 


N v’ 

’*su v 


k 


,v 


N c 


v 


M 

C 


(30) 
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*0 - N s*v m (31) 

where 

U o * {^o » ^0« ^0 1 anc ^ ® 0 “ l^x” * ®y" i 

v" T - {^oi * ^01 * ^01 * ®x"1 * ®y"1 ■ ®z"1 * ^03 > ^03- ^03 » ®x"3 * ^y"3* ®z"3) 

V ’ T * {^04 » — • ^09 • ^04* — » ^09 « ^04* ^05 » ^06 1 
v "c " i^Oc • ^0c) 

The matrix of interpolation functions, N su , is dependent on yV^ s - and (i - 1, 2, 3). 
The shape function matrix, fsJ su , involves ^ si (i - 4 f 9) and (i * 4, .... 6). The 
matrix N sc contains only yV^sc. The shape matrix, N s0f is composed of the area 
coordinates |j (i - 1, 3). Although 0 2 « k (k - 1, 2, 3) do not exist, they are included in the 
formulation because the transformation to body-fixed coordinates leads to non-zero 
rotations in all axes. 

The intra-edge displacements included in v" are expressed in terms of corner 
displacements, v", through an appropriate transformation 

v" = L v" (32) 

The transformation matrix, L, is obtained by imposing constant strain conditions along 
the three edges of the element. The explicit form of these constraints and the detailed 
form of L (in partitioned form) are given by Tessler [ 1 1 ]. 

Substituting from Eq. (4) into Eq. (30) and combining the resulting equation with Eq. 
(31) lead to 



"where N su - N su + i^ su L. This relationship can be taken as the basis in formulating the 
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strain-displacement transformation matrix, £3 S , in the form 

- (^ SU , f3 s g] 


( 34 ) 


where 68^ and <J ? s9 are defined such that the resultant strain vector, t e L , given by Eq. 
(26) can be obtained from 


t°L “ [®su- 


r ^ 

u o 

®o 


(35) 


This form of vector t e L allows the first integral in Eq. (24) to be rewritten as 





T 





■ 

▼ 

5v" 



k L11 

k L12 


v" ' 

5 x e(_ C t e L d - 4 



> 



4 

> 


. 5V '= 

- 


k T 

*112 

k L22 


V'c 


(36) 


with 


'Ll 1 


l L12 


(n^Lc ® su n; u + n;^j u c ® s9 n s9 

«A 

+ fi8 su N; u + Nj fl #J 9 C <3 s9 N s9 )d ‘a 

( N su®I u C & su N sc + Nj 0 <Zj 9 C & su N sc )d 'a 
l A 


and 


'L22 




Nj c Sj u C fi? su N sc d <a 


Assuming external or internal loads are not applied at the center node of the element, 
the stiffness matrix, k L , can be condensed statically to a smaller sized matrix, k L , 


k L = 


'•L11 


- k |_12 k L22 k L12 


(37) 


where k L represents the condensed stiffness matrix. It consists of six degrees of 
freedom, including the 6 Z artificial degree of freedom at each corner node. 
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The second integral in Eq. (24) can be rewritten in terms of the total nodal vector 
as 

5 t^NL * s m d ' a * * v"T k g v" (38) 

Jt A 

where the geometric stiffness matrix, k g , is defined by 


'S m £8 N, 


mu 


da 


(39) 


m which the matrix differential operator £8 and ‘S,,, are explicitly given by Madenci 
and Barut [ 1 3] . The tangential stiffness matrix, k T , is defined by 


k T = k i. + k g (40) 

By substituting from Eqs. (25) and (35), the third integral in Eq. (24) can be 
rewritten in the form 


5 t eL l s d'a - 5 v" T k[ v " (41) 

Jt A 

The vector of nodal displacements, v"*, corresponds to local deformations between the 
initial and the updated configurations of the elements. It is determined by subtracting 
the rigid-body motion from the total displacements. The deformational out-of-plane 
nodal rotations are determined by the methods provided by Rankin and Brogan [14], 
With the use of Eqs. (36), (38), and (41), in conjunction with Eq. (11), the second 
integral, representing the virtual work by internal stresses in the kth element, is 


expressed in matrix form as 
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where 

k' T - Z'(k[ + k fl )Z' T and k[_ - Z’k L Z' T 

Virtual work by external forces. The virtual work due to nodal forces and 
moments at the ith node of the kth element is expressed as 

hW y - SR^Fj + (Sw T + ^JMj (43) 

in which F; and M; are the nodal force and moment vectors, respectively, at the ith 
node. The position of the ith node in the kth element with respect to the origin of the 
inertial frame is specified by Rj. The virtual rigid-body rotation of this element is 
represented by ha?-, and its elastic rotation at the ith node by a vector, 50,. Substituting 
for 5Rj from Eq. (5) into Eq. (43) and using the property of the transformation matrix, T, 
result in 

hW s = 6RjF| + 5«/ r T(r' j F' j + M',) + 5ufF’j + M'J M', (44) 

Adding the contribution of each node and utilizing the definition of v'^ result in the 
following expression for virtual work due to external loads in the kth element: 



where 

3 3 

F - Fi 1 M- T XViF’ +M’,) 
i=1 i-1 


and 
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P ,J -{P'J P'l P'J} . With P’i-iFj Mft 

Finite element equations of motion. Substituting Eqs. (21), (42), and (45) into the 
expression for virtual work (Eq. 1) and requiring the virtual quantities 6R 0 , 5 w, and W 
to be arbitrary result in the coupled and highly nonlinear equations of motion between 
time t and t + At for the finite element: 


M(b)b + f (b) - f(b) 


(46) 


The solution vector, b, the mass matrix, M, and the vector, f, arising from external and 
gyroscopic forces are defined as 





” 


r > 


Ro 


M rr Mfw M Rv . 


p - ft - h 

b « < 

co- 

►; M - 

^ u> w- M^ v . 

; f - • 

M"- 9t - 92 


v' 


Sym. M vV 


P’ - h, - h 2 


The vector due to internal stresses, f ff , is given by 

f = K T (b)b + 'f a (47) 

where 



0 0 0 


0 

K t - 

0 0 0 

and t f or = < 

0 


1 

_h- 

o 

O 

i 


k^V" 


The vector l f a is known from the previous increment between time t - At and t. 

The equations of motion (46) for each element are assembled by usual techniques 
to form the global equations of motion. Their solution provides the motion of the origin 
of the body-fixed frame with respect to the inertial frame and the elastic deformations 
of the body with respect to the body-fixed frame. 
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Numerical Results 

Solution of the nonlinear global equations of motion is achieved by means of the 
Newton-Raphson method in conjunction with the trapezoidal time integration scheme. 
The construction of the solution vector b at time t + At involves iterations that begin 
with the known solution vector ’b at time t. The ith iteration of the solution vector b ( is 
decomposed as 

bj - b j _ 1 + Ab (48) 

In constructing the solution vector b, the ith iteration of the solution method puts the 
equations of motion (46) in the form 

M(bj )t>j + f ff (bj) - f(bj) (49) 


In order to gain computational efficiency in the construction of the mass and geometric 
stiffness matrices, M and k g , respectively, the shape matrix, N m , is a reduced form of 
the shape matrix for the shallow shell element. The explicit form of N m can be found in 
Ref. [ 1 5] . The validity of the present analysis is established by considering a rotating 
beam subjected to a spin-up maneuver [12]. In this study, the beam is modeled as a 
narrow plate with triangular shallow shell elements with the following boundary 
conditions: 

* 0 (o. y'. o) ■ y'» “ ^o(®- y'' o) ** o 


dx' 


(0, y\ 0) - 


dy' 


(0, y', 0) - 0 


The spin-up maneuver for this problem is specified as 
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where J2 0 - rr/ 1 0 rad/sec and T 0 - 1 sec. The Young's and shear moduli for the beam 
are given as E - 1.44 x 10® lb/ ft 2 and G - 5.54 x 10 7 lb/ft 2 . The cross-sectional area 
and the length of the beam are A - 1 ft 2 and L - 100 ft, respectively. The moment of 
inertia about the planar axes is specified as I - 0.08333 ft 4 . The beam has a mass 
density of p - 5.22 slug/ft 3 . Since the beam is modeled as a plate, its thickness and 
width are computed as h - 1 ft and W - 1 ft based on the values given for the moment 
of inertia and the cross-sectional area for the beam. A comparison between the 
considerably large transverse and axial displacements at the tip of the plate obtained 
from the present analysis and those presented by Christensen and Lee is given in Fig- 
3. The sequence of the motion with respect to the inertial coordinate system is 
depicted in Fig. 4. 

The capability of the present analysis is further demonstrated by considering a 
laminated, cantilever blade tilted at 45° and 90° angles with respect to the (x-y) plane. 
The planform dimensions of the blade are specified by L - 2.5 m and W - 0.2 m. The 
blade thickness of h - 0.04 m results from a laminate lay-up of [ 0° /905 ] s . Each lamina 
has the material properties E u - 181 GPa, E T - 10.3 GPa, G LT - 7.17 GPA, and i> LT = 
0.28, with mass density p - 1600 kg/m 3 . The applied spin-up maneuver is the same as 
that specified for the previous problem, with Q 0 - 10 jt rad/sec and T 0 - 0.05 sec. The 
influence of tilt angle on the tip displacements is illustrated in Fig. 5. As shown in this 
figure, the difference in maximum tip deflections is rather significant. When the blade is 
tilted at a 45° angle, the difference between the transverse deflections at the left and 
right corners of the blade indicates that the tip of the blade undergoes twisting. These 
corner displacements are presented in Fig. 6. With respect to the inertial frame, the 
sequence of deformations for a blade tilted at angles of 0° and 45° is illustrated in Fig. 


7. 
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Conclusions 

A variational formulation that accounts for changes in motion-induced stiffness 
while determining the structural response of moderately thick laminated shells within 
the realm of a multibody dynamics formalism has been presented. This methodology 
also provides the inertial forces as part of the solution for corresponding external 
forces. Although demonstrated herein on simple geometric configurations, this analysis 
is general enough to treat curved laminates. This formulation represents an 
advancement in state-of-the-art analysis of rigid-body dynamics coupled with nonlinear 
elastic deformations. 
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FIGURE CAPTIONS 

Figure 1. Kinematics of a point in an element due to rigid-body motion and local 
deformation. 

Figure 2. Element coordinate system and degrees of freedom for a triangular shallow 
shell element. 

Figure 3. Transverse and axial tip displacements along the center line of a rotating 
narrow plate 

Figure 4. A narrow plate undergoing a sequence of deformations in reference to an 
inertial frame. 

Figure 5. Transverse tip displacements along the center line of a rotating laminated 
blade. 


Figure 6. Transverse tip displacements at the corners of a rotating laminated blade. 

Figure 7. Sequence of deformations for a blade tilted at angles of 90° and 45°. 
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